source('00_util_scripts/mod_bplot.R')
library(Peptides)

pept <- read_tsv('virtual_screen_v3/10aa_trpv3-fpp.tsv')

g1 <- pept |>
  mutate(pi = pI(aaseq),
         hydrophobicity = hydrophobicity(aaseq)) |>
  ggplot(aes(x = ID, y = hydrophobicity)) +
  geom_path() +
  geom_point() +
  theme_bw() +
  geom_hline(yintercept = 6.5, linetype = 'dashed') +
  labs(title = 'Electric charges of candidate TRPV3-binding peptides',
       x = 'rank', y = 'PI')

g1

hiv_tat <- 'RRRQRRKKRG'

pI(hiv_tat)

hiv_tat_natural <- str_split_1(hiv_tat, '') |> rev() |> str_c(collapse = '')

pI(hiv_tat_natural)

r7 <- strrep('R', 7)

k7 <- strrep('K', 7)

penetratin <- 'RQIKWFQNRRMKWKK'
ptd4 <- 'YARAAARQARA'
pI(penetratin)

cpps <-
tibble(name = c('TAT(48-57)', 'PTD4', 'R7', 'K7'),
       aaseq = c(hiv_tat, ptd4, r7, k7),
       pi = pI(aaseq),
       ID = c(5,30,55,80))

g1 + geom_point(data = cpps) +
  geom_label_repel(data = cpps, aes(label = name))

pept |>
  mutate(tat_aa = str_c(hiv_tat_natural, aaseq),
         aa_tat = str_c(aaseq, hiv_tat_natural),
         ptd4_aa = str_c(ptd4, aaseq),
         aa_ptd4 = str_c(aaseq, ptd4),
         r7_aa = str_c(r7, aaseq),
         aa_r7 = str_c(aaseq, r7)) |>
  pivot_longer(-1, names_to = 'type', values_to = 'aaseq') |>
  write_csv('virtual_screen_v3/cpp_10aa_trpv3-fpp.csv')

pept <- read_csv('virtual_screen_v3/cpp_10aa_trpv3-fpp.csv')

pept |>
  filter(type == 'tat_aa') |>
  mutate(smiles = aaSMILES(aaseq)) |>
  pull(smiles) |>
  head(1)

# prepare protenix inputs ---------
library(jsonlite)
prtnx <- read_json('virtual_screen_v3/inputs.json', simplifyVector = T)

pept |>
  filter(type == 'tat_aa') |>
  pull(aaseq, name = ID) |>
  iwalk(\(x, idx){
    prtnx$sequences[[1]]$proteinChain$sequence[[2]] <- x
    prtnx |> write_json(str_glue('virtual_screen_v3/results/inputs_{idx}.json'))
    })

tat_pept <- pept |>
  filter(type == 'tat_aa')

prtnx_prot_seq <-
tat_pept$aaseq |>
  map(\(x)list(proteinChain = list(count = 1, sequence = x, modification = list())))

prtnx[[1]]$sequences <- c(prtnx[[1]]$sequences, prtnx_prot_seq[-1])

prtnx |> glimpse()

prtnx |> write_json('virtual_screen_v3/msa_input.json', auto_unbox = T)

# prepare Alphafold3 input ---------
allrank <- read_csv('virtual_screen_v3/results/model3_all_rank.csv')

dscf_300 <- allrank |>
  filter(!is.na(aaseq))

af3_b1 <- dscf_300 |>
  slice_min(DSCF_rank, n = 30) |>
  pull(aaseq, name = ID) |>
  imap(\(x, idx){
    prtnx$name <- idx
    prtnx$use_esm <- NULL
    prtnx$use_msa <- NULL
    prtnx$covalent_bonds <- NULL
    prtnx$constraint <- NULL
    prtnx$sequences[[1]]$proteinChain$sequence[[2]] <- x
    prtnx
  })

af3_b1 |> glimpse()

af3_b1 <- unname(af3_b1)

af3_b1 |> list_simplify() |>
  write_json('virtual_screen_v3/results/af3_input_b1.json',
         auto_unbox = T, pretty = T)

dscf_300 <- dscf_300 |>
  mutate(DSCF_rank = min_rank(-DSCF_score))

betw2 <- seq.int(30,300,30)
betw1 <- betw2 - 29

list(betw1, betw2) |>
  pmap(\(x,y){dscf_300 |>
      filter(between(DSCF_rank, left = x, right = y)) |>
      pull(aaseq, name = ID) |>
      imap(\(x, idx){
        prtnx$name <- idx
        prtnx$use_esm <- NULL
        prtnx$use_msa <- NULL
        prtnx$covalent_bonds <- NULL
        prtnx$constraint <- NULL
        prtnx$sequences[[1]]$proteinChain$sequence[[2]] <- x
        prtnx
      }) |>
      unname() |>
      list_simplify()}) |>
  iwalk(\(x, idx){
    write_json(x, str_glue('virtual_screen_v3/results/af3_input_b{idx}.json'),
               auto_unbox = T, pretty = T)})

# AF3 output ---------
af3_b1o <- list.files('~/append-ssd/ringtail-test/data/af3_predict/',
           '.out', recursive = T, full.names = T) |>
  set_names() |>
  map(read_file)

af3_b1o_score <- af3_b1o |>
  list_simplify() |>
  str_extract(' -.+\\d') |>
  str_squish() |>
  as.double()

names(af3_b1o) |> str_subset('peptide', negate = T)

af3_rank <-
tibble(file = names(af3_b1o) |> basename(),
       score = af3_b1o_score,
       id = str_extract(file, '(?<=fold_).+(?=_model)') |> str_remove('^trpv3_') |>
         str_to_upper() |> str_replace('PEPTIDE', 'Peptide'),
       conf_rank = str_extract(file, '\\d(?=.cif)'),
       model = case_match(conf_rank, '0' ~ 'AF3',
                          '1' ~ 'Protenix',
                          '2' ~ 'Chai1',
                          .default = NA)) |>
  na.omit() |>
  mutate(rank = row_number(score), .by = model)

af3_rank |>
  ggplot(aes(rank, -score)) +
  geom_point() +
  facet_wrap(~model) +
  theme_bw() +
  scale_x_reverse()

af3_rank |>
  write_csv('virtual_screen_v3/results/alphafold3_300cpp-peptide_TRPV3_score.csv')

af3sv_rank <- af3_rank |>
  filter(conf_rank == '0')

af3_top <- af3sv_rank |>
  slice_min(score, n=8)

af3sv_rank |>
  ggplot(aes(rank, -score)) +
  geom_point() +
  geom_label_repel(data = af3_top, aes(label = id), nudge_x = -20, nudge_y = 1) +
  geom_point(data = af3_top, color = 'red') +
  labs(y = 'Binding energy (-Kcal/mol)', x = 'Rank',
       subtitle = 'Abramson et al., 2024 (Google)',
       title = 'AlphaFold3 predicted binding score with TRPV3') +
  scale_x_reverse() +
  theme_bw()

prtnx_rank <- af3_rank |>
  filter(conf_rank == '1')

prtnx_top <- prtnx_rank |>
  slice_min(score, n=8)

prtnx_rank |>
  ggplot(aes(rank, -score)) +
  geom_point() +
  geom_label_repel(data = prtnx_top, aes(label = id), nudge_x = -20, nudge_y = 1) +
  geom_point(data = prtnx_top, color = 'red') +
  labs(y = 'Binding energy (-Kcal/mol)', x = 'Rank',
       subtitle = 'Chen et al., 2025 (ByteDance)',
       title = 'Protenix predicted binding score with TRPV3') +
  scale_x_reverse() +
  theme_bw()

chai_rank <- af3_rank |>
  filter(conf_rank == '2')

chai_top <- chai_rank |>
  slice_min(score, n=8)

chai_rank |>
  ggplot(aes(rank, -score)) +
  geom_point() +
  geom_label_repel(data = chai_top, aes(label = id), nudge_x = -20, nudge_y = 1) +
  geom_point(data = chai_top, color = 'red') +
  labs(y = 'Binding energy (-Kcal/mol)', x = 'Rank',
       subtitle = 'Boitreaud et al., 2024 (Chai Discovery)',
       title = 'Chai1 predicted binding score with TRPV3') +
  scale_x_reverse() +
  theme_bw()

## robust? ----------
af3_wide_rank <- af3_rank |>
  pivot_wider(id_cols = id, names_from = model,
              values_from = c(rank, score))

af3_wide_rank |>
  rowwise() |>
  mutate(max_rank = max(rank_AF3, rank_Protenix, rank_Chai1)) |>
  arrange(max_rank) |>
  head(97) |>
  ungroup() |>
  summarise(AlphaFold3 = sum(rank_AF3 <= 96),
            Protenix = sum(rank_Protenix <= 96),
            Chai1 = sum(rank_Chai1 <= 96)) |>
  pivot_longer(1:3) |>
  ggplot(aes(name, value, fill = name)) +
  geom_col() +
  geom_text(aes(label = value), nudge_y = -5, color = 'white') +
  expand_limits(y = 96) +
  labs(title = 'Peptides belong to top-96 of single-model',
       x = 'Model', y = 'Count', fill = 'Model') +
  theme_bw(base_size = 14)

intersect_n <- function(x){
  af3_wide_rank |>
    filter(rank_AF3 < x, rank_Protenix < x, rank_Chai1 < x,
           str_detect(id, 'PTD4|TAT|R7')) |>
    nrow()
}

allrank |>
  select(ID, aaseq) |>
  right_join(af3_wide_rank, join_by(ID == id)) |>
  write_csv('virtual_screen_v3/results/af3_based_model3_wide_rank.csv')

intsct_af33 <- 1:300 |>
  map_int(intersect_n, .progress = T)

tibble(threshold = 1:300, size = intsct_af33) |>
  ggplot(aes(threshold, size)) +
  geom_point(aes(color = size <= 96)) +
  annotate(geom = 'segment', x = 169, y = 0, xend = 169, yend = 96, linetype = 'longdash') +
  annotate(geom = 'segment', x = 0, y = 96, xend = 169, yend = 96, linetype = 'longdash') +
  theme_bw() +
  labs(x = 'Rank threshold', y = 'Intersection size', color = 'Top 96',
       title = 'Intersection of top rank peptides from 3 AF3-based models')

# prepare DeepSCFold input ---------
library(Biostrings)

tat_pept |>
  mutate(ID = str_c(type, ID)) |>
  pull(aaseq, name = ID) |>
  AAStringSet() |>
  writeXStringSet('tat_aa96.fasta')

pept |>
  mutate(ID = str_c(type, ID, sep = '_')) |>
  pull(aaseq, name = ID) |>
  AAStringSet() |>
  writeXStringSet('peptide762.fasta')

allq <-
readAAStringSet('~/append-ssd/project_dev/deepscfold/ppi_pred/all_v3_query762.fasta')

allq

## PTD4-1aa effect on binding ----------
widerank <- read_csv('virtual_screen_v3/results/af3_based_model3_wide_rank.csv')

widerank |>
  filter(str_detect(ID, 'PTD4')) |>
  mutate(aaseq_n1 = str_replace(aaseq, 'YARAAAR', 'ARAAAR'),
         aaseq_c1 = str_replace(aaseq, 'ARQARA', 'ARQAR')) |>
  select(c(ID, contains('aaseq'))) |>
  pivot_longer(-1) |>
  mutate(ID = str_c(ID, name)) |>
  pull(value, name = ID) |>
  AAStringSet() |>
  writeXStringSet('ptd4_minus1.fasta')

ptd4_m1 <-
  read_csv('~/append-ssd/project_dev/deepscfold/ppi_pred/ptd4_minus1_out/final_scores.csv')

ptd4_m1 |>
  mutate(ID = str_extract(QuerySeqID, '.+aaseq'),
         type = str_extract(QuerySeqID, 'seq.*')) |>
  pivot_wider(names_from = type, values_from = Score, id_cols = ID) |>
  mutate(delta_c1 = seq_c1 - seq, delta_n1 = seq_n1 - seq,
         fusion = ifelse(str_detect(ID, '^PTD4'), 'N-PTD4', 'PTD4-C')) |>
  ggplot() +
  geom_boxplot(aes(x = 'Delta_N', y = delta_n1)) +
  geom_boxplot(aes(x = 'Delta_C', y = delta_c1)) +
  facet_wrap(~fusion)

ptd4_m1 <-
read_csv('~/append-ssd/project_dev/deepscfold/ppi_pred/ptd4_minus1_out_rev/final_scores.csv')

ptd4_m1 |>
  mutate(ID = str_extract(TargetSeqID, '.+aaseq'),
         type = str_extract(TargetSeqID, 'seq.*')) |>
  pivot_wider(names_from = type, values_from = Score, id_cols = ID) |>
  mutate(delta_c1 = seq_c1 - seq, delta_n1 = seq_n1 - seq,
         fusion = ifelse(str_detect(ID, '^PTD4'), 'N-PTD4', 'PTD4-C')) |>
  ggplot() +
  geom_boxplot(aes(x = 'Delta_N', y = delta_n1)) +
  geom_boxplot(aes(x = 'Delta_C', y = delta_c1)) +
  facet_wrap(~fusion)

## deepscfold result -------
dsf_path <- '~/append-ssd/project_dev/deepscfold/ppi_pred/'

res_model1 <-
read_csv(str_glue('{dsf_path}output_762_trpv3_model1/final_scores.csv'))

res_model1 <- res_model1 |>
  mutate(rank = rank(Score),
         QuerySeqID = ifelse(str_detect(QuerySeqID, 'TRPV3'), 'TRPV3',
                             QuerySeqID) |>
           str_to_upper() |> str_replace('AA', 'Peptide') |> str_remove('SEQ'))

top6_model1 <- res_model1 |>
  slice_max(Score, n = 8)

known_model1 <- res_model1 |>
  filter(!str_detect(QuerySeqID, 'Peptide'))

res_model1 |>
  ggplot(aes(rank, Score)) +
  geom_point() +
  scale_y_log10() +
  geom_point(data = top6_model1, color = 'red') +
  geom_label_repel(data = top6_model1, aes(label = QuerySeqID), nudge_x = -100) +
  theme_bw() +
  labs(title = 'DeepSCFold predicted binding score with TRPV3',
       subtitle = 'Hou et al., 2025 (Westlake Univ)')

res_model2 <-
  read_csv(str_glue('{dsf_path}output_762_trpv3_model2/final_scores.csv'))

res_model2 <- res_model2 |>
  mutate(rank = rank(Score),
         QuerySeqID = ifelse(str_detect(QuerySeqID, 'TRPV3'), 'TRPV3',
                             QuerySeqID) |>
           str_to_upper() |> str_replace('AA', 'Peptide') |> str_remove('SEQ'))

top_model2 <- res_model2 |>
  slice_max(Score, n = 8)

known_model2 <- res_model2 |>
  filter(!str_detect(QuerySeqID, 'Peptide'))

res_model2 |>
  ggplot(aes(rank, Score*100)) +
  geom_point() +
  #scale_y_log10() +
  geom_point(data = top_model2, color = 'red') +
  geom_label_repel(data = top_model2, aes(label = QuerySeqID), nudge_x = -100) +
  theme_bw() +
  labs(title = 'Autodock-Crankpep predicted binding score with TRPV3',
       subtitle = 'Zhang et al., 2019 (Scripps Institute)',
       y = '-Free energy (kCal/mol)')

# prepare PIPR input -----------
allq <-
  readAAStringSet('~/append-ssd/project_dev/deepscfold/ppi_pred/all_v3_query762.fasta')

allq_char <- allq |> as.character()

tibble(p1 = allq_char, p2 = allq_char[['human_TRPV3']]) |>
  write_tsv('pipr_v3_678.tsv', col_names = F)

## PIPR result -----
pipr678 <-
read_delim('~/append-ssd/project_dev/seq_ppi/binary/model_sun_human/lasagna/v3_678_epoch100_pred.tsv')

pipr678$ID <- names(allq_char) |> str_remove('human_') |>
  str_to_upper() |> str_replace('AA', 'Peptide')

pipr678$Prediction |> table()

pipr678_rank <- pipr678 |>
  mutate(rank = rank(Probability))

pipr678_top <- pipr678_rank |>
  slice_max(Probability, n = 8, with_ties = F)

pipr678_rank |>
  ggplot(aes(rank, Probability)) +
  geom_point() +
  scale_y_log10() +
  geom_point(data = pipr678_top, color = 'red') +
  geom_label_repel(data = pipr678_top, aes(label = ID), nudge_x = -100) +
  labs(title = 'PIPR predicted binding score with TRPV3',
       subtitle = 'Chen et al., 2019 (UCLA)') +
  theme_bw()

## PIPR RAM --------
pipr_gpu |>
  mutate(timestamp = as_datetime(timestamp),
         memory_use = str_extract(`utilization.memory [%]`, '\\d+') |> as.double(),
         memory_allocate = str_extract(`memory.used [MiB]`, '\\d+') |> as.double()) |>
  filter(timestamp < as_datetime('2025-6-24 21:50:00')) |>
  ggplot(aes(timestamp)) +
  #geom_path(aes(y = memory_use)) +
  geom_path(aes(y = memory_allocate)) +
  facet_grid(rows = vars(index))

# prepare boltz2 input -----------
library(yaml)

allrank <- read_csv('virtual_screen_v3/results/model3_all_rank.csv')

pipr_inp <- read_tsv('virtual_screen_v3/pipr_v3_678.tsv', col_names = F)

full_trpv3 <- pipr_inp$X2[[1]]

nchar(full_trpv3)

ptd4_pept_30 <- allrank |>
  na.omit() |>
  slice_max(DSCF_score) |>
  pull(aaseq)

list(sequences = list(list(protein = list(id = 'A',
                                     sequence = full_trpv3)),
                      list(protein = list(id = 'B',
                                     sequence = ptd4_pept_30)))) |>
  write_yaml('virtual_screen_v3/results/ptd4_pept_30.yaml')

write_boltz_yaml <- function(peptide, id, target_seq, yaml_dir, msa_dir, verbose = F){
  require(yaml)
  dir.create(yaml_dir, showWarnings = F, mode = '0776')
  dir.create(msa_dir, showWarnings = F, mode = '0776')
  yaml_path <- str_glue('{yaml_dir}/{id}.yaml')
  msa_path <- str_glue('{msa_dir}/{id}.csv')
  list(sequences = list(list(protein = list(id = 'A',
                                            msa = str_glue('./{msa_dir}/target_prot_msa.csv'),
                                            sequence = target_seq)),
                        list(protein = list(id = 'B',
                                            msa = msa_path,
                                            sequence = peptide)))) |>
    write_yaml(yaml_path)
  if (verbose) message(str_glue('Yaml file written to {yaml_path}.'))
  tibble(key = c(0,-1), sequence = peptide) |>
    write_csv(msa_path)
  if (verbose) message(str_glue('MSA file written to {msa_path}.'))
  }

allrank |>
  na.omit() |>
  slice_max(DSCF_score, n = 300) |>
  mutate(peptide = aaseq, id = ID,
         yaml_dir = '~/append-ssd/learn/boltz/v3_pept_yaml',
         .keep = 'none') |>
  pwalk(write_boltz_yaml)

## boltz2 output ----------
boltz_json <-
'~/append-ssd/learn/boltz/v3_pept_msa_out_dummy/boltz_results_v3_pept_yaml/predictions/' |>
  list.files(pattern = 'json', recursive = T, full.names = T)

boltz_conf <- boltz_json |>
  map(\(x)read_json(x, simplifyVector = T) |> as_tibble() |> head(1)) |>
  list_c()
  
boltz_conf |>
  select(confidence_score, ptm, iptm, complex_plddt, complex_iplddt) |>
  pivot_longer(everything()) |>
  ggplot(aes(x = 0, y = value)) +
  geom_boxplot() +
  geom_jitter(height = 0, width = .1) +
  facet_wrap(~name) +
  theme_bw() +
  labs(title = 'Boltz2 prediction QC of pocket dummy settings (n=300)',
       x = 'QC score type') +
  theme(axis.text.x = element_blank())

boltz_out <-
  '~/append-ssd/learn/boltz/v3_pept_msa_out/boltz_results_v3_pept_yaml/' |>
  list.files(pattern = 'prodigy.out$', recursive = T, full.names = T) |>
  read_lines()

boltz_parsed <-
tibble(out = boltz_out, type = ifelse(str_ends(out, '\\)'), 'id', 'score')) |>
  mutate(parsed = str_extract(out, '[:graph:]+model|[:graph:]+\\d$'))

boltz_id <- boltz_parsed |>
  filter(type == 'id') |>
  mutate(id = str_remove(parsed, '_model'), .keep = 'none')

boltz_score <- boltz_parsed |>
  filter(type == 'score') |>
  mutate(score = as.double(parsed), .keep = 'none') |>
  bind_cols(boltz_id) |>
  mutate(rank = row_number(score))

boltz_top <- boltz_score |>
  slice_min(score, n = 8)

boltz_score |>
  ggplot(aes(rank, -score, label = id)) +
  geom_point() +
  scale_x_reverse() +
  theme_bw() +
  geom_label_repel(data = boltz_top, aes(label = id), nudge_x = -20, nudge_y = 1) +
  geom_point(data = boltz_top, color = 'red') +
  labs(y = 'Binding energy (-kCal/mol)', subtitle = 'Passaro et al., 2025 (MIT)',
       title = 'Boltz2 pocket predicted binding score')

# pepmimic pept into boltz2 --------------
pep_raw <- list.files('~/append-ssd/data_lfs/results/', full.names = T) |>
  read_delim(col_names = pepcol, id = 'file') |>
  mutate(metric = str_extract(file, '(?<=results//).+(?=.txt)'), .keep = 'unused')

pep_raw <- pep_raw |>
  pivot_wider(names_from = metric, values_from = score)

pep_filter <- pep_raw |>
  filter(foldx_dG < -5, pyrosetta_dG < -10, interface_hit > 2) |>
  mutate(rank_foldx = rank(foldx_dG), rank_pyrosetta = rank(pyrosetta_dG),
         rank_mean = rank_foldx + rank_pyrosetta) |>
  slice_min(rank_mean, n = 500)

cd19_ecm <- 'PEEPLVVKVEEGDNAVLQCLKGTSDGPTQQLTWSRESPLKPFLKLSLGLPGLGIHMRPLAIWLFIFNVSQQMGGFYLCQPGPPSEKAWQPGWTVNVEGSGELFRWNVSDLGGLGCGLKNRSSEGPSSPSGKLMSPKLYVWAKDRPEIWEGEPPCLPPRDSLNQSLSQDLTMAPGSTLWLSCGVPPDSVSRGPLSWTHVHPKGPKSLLSLELKDDRPARDMWVMETGLLLPRATAQDAGKYYCHRGNLTMSFHLEIT'

pep_filter$seq |>
  walk2(pep_filter$id, \(x,y)write_boltz_yaml(
    peptide = x, id = y, target_seq = cd19_ecm,
    yaml_dir = '~/append-ssd/learn/boltz/cd19_pepmimic_boltz2_yaml',
    msa_dir = 'cd19_pepmimic_boltz2_msa', verbose = T)
  )

## boltz2 output ----------
boltz_json <-
  '~/append-ssd/learn/boltz/cd19_pepmimic_ecd_5recycle_all/boltz_results_cd19_pepmimic_boltz2_yaml/predictions/' |>
  list.files(pattern = 'json', recursive = T, full.names = T)

boltz_conf <- boltz_json |>
  map(\(x)read_json(x, simplifyVector = T) |> as_tibble() |> head(1)) |>
  list_c()

boltz_conf |>
  select(confidence_score, ptm, iptm, complex_plddt, complex_iplddt) |>
  pivot_longer(everything()) |>
  ggplot(aes(x = 0, y = value)) +
  geom_boxplot() +
  geom_jitter(height = 0, width = .1) +
  facet_wrap(~name, scales = 'free_y') +
  theme_bw() +
  labs(title = 'Boltz2 prediction QC of CD19 pepmimic peptides (n=500)',
       x = 'QC score type') +
  theme(axis.text.x = element_blank())

# reproducibity among models -------
g1 <- res_model1 |>
  left_join(res_model2, join_by(QuerySeqID)) |>
  ggplot(aes(rank.x, rank.y)) +
  geom_point() +
  stat_cor() +
  geom_smooth(method = 'lm') +
  theme_bw() +
  labs(x = 'Rank by DeepSCFold', y = 'Rank by Autodock-Crankpep',
       title = 'Correlation of result ranking')

g2 <- res_model1 |>
  inner_join(pipr678_rank, join_by(QuerySeqID == ID)) |>
  ggplot(aes(rank.x, rank.y)) +
  geom_point() +
  stat_cor(geom = 'label', alpha = .5) +
  geom_smooth(method = 'lm') +
  theme_bw() +
  labs(x = 'Rank by DeepSCFold', y = 'Rank by PIPR')

g3 <- res_model2 |>
  inner_join(pipr678_rank, join_by(QuerySeqID == ID)) |>
  ggplot(aes(rank.x, rank.y)) +
  geom_point() +
  stat_cor(geom = 'label', alpha = .5) +
  geom_smooth(method = 'lm') +
  theme_bw() +
  labs(x = 'Rank by DeepSCFold', y = 'Rank by PIPR',
       title = 'Correlation of result ranking')

g1 + g2 + g3

z21_normalize <- function(x) {
  scale(x, center = min(x), scale = max(x))[,1]}

res_model1 |>
  mutate(scaled_score = z21_normalize(Score)) |>
  ggplot(aes(rank, scaled_score)) + geom_point()

pipr_scale <- pipr678_rank |>
  mutate(ID, PIPR_score = z21_normalize(Probability), .keep = 'none')

dsf_scale <- res_model1 |>
  mutate(ID = QuerySeqID, DSCF_score = z21_normalize(Score), .keep = 'none')

dsf2_scale <- res_model2 |>
  mutate(ID = QuerySeqID, ADCP_score = z21_normalize(Score), .keep = 'none')

merge_scale <- pipr_scale |>
  mutate(ID = str_remove(ID, 'SEQ')) |>
  inner_join(dsf_scale) |>
  inner_join(dsf2_scale) |>
  pivot_longer(-ID) |>
  summarise(value = mean(value), .by = ID) |>
  mutate(rank = rank(value))

merge_top <- merge_scale |>
  slice_max(value, n = 8)

merge_red <- merge_scale |>
  slice_max(value, n = 101)

merge_scale |>
  ggplot(aes(678-rank, value)) +
  geom_point() +
  geom_vline(xintercept = 101, linetype = 'dashed') +
  scale_x_reverse() +
  scale_y_log10() +
  geom_label_repel(data = merge_top, aes(label = ID),
                   nudge_x = -10, nudge_y = .16) +
  geom_point(data = merge_red, color = 'red') +
  labs(x = 'Rank', y = 'Binding score',
       title = 'Average predicted binding score with TRPV3') +
  theme_bw()

top96_model3 <- pept |>
  mutate(ID = str_to_upper(type) |> str_replace('AA', 'Peptide') |>
           str_c('_', ID) |> str_remove('SEQ'),
         type = str_remove(ID, '_\\d+$')) |>
  right_join(merge_scale) |>
  mutate(rank = rank(-value)) |>
  arrange(desc(value)) |>
  write_csv('virtual_screen_v3/results/model3_scored_peptide678.csv')

top96_model3 |>
  filter(!is.na(type)) |>
  ggplot(aes(675-rank, color = type)) +
  geom_density() +
  theme_bw() +
  scale_color_viridis_d(option = 'turbo', begin = .1)

top96_model3 |>
  filter(!is.na(type)) |>
  mutate(type = str_to_upper(type) |> str_replace('AA', 'Peptide') |>
          str_remove('SEQ')) |>
  dplyr::count(type) |>
  mutate(type = fct_reorder(type, n)) |>
  ggplot(aes(type, n, fill = type)) +
  geom_col() +
  theme_bw() +
  scale_color_viridis_d(option = 'turbo', begin = .1) +
  labs(title = 'CPP fusion modes in top 96 scored peptides',
       y = 'Count')

## intersect of tops -----
inter_scales <- pipr_scale |>
  mutate(ID = str_remove(ID, 'SEQ')) |>
  inner_join(dsf_scale) |>
  inner_join(dsf2_scale)

inter_ranks <- inter_scales |>
  mutate(mean_score = mean(c(PIPR_score, DSCF_score, ADCP_score)), .by = ID) |>
  mutate(PIPR_rank = min_rank(-PIPR_score),
         DSCF_rank = min_rank(-DSCF_score),
         ADCP_rank = min_rank(-ADCP_score),
         mean_rank = min_rank(-mean_score))

pept |>
  mutate(ID = str_to_upper(type) |> str_replace('AA', 'Peptide') |>
           str_c('_', ID) |> str_remove('SEQ'),
         type = str_remove(ID, '_\\d+$')) |>
  right_join(inter_ranks) |>
  write_csv('virtual_screen_v3/results/model3_all_rank.csv')

# A431 in vitro cytokine ------------
a431_cytk <- read_csv('virtual_screen_v3/data/peptide96-cytokine.csv')

high_cytk <- a431_cytk |>
  mutate(rank_il6 = percent_rank(IL6), rank_ccl20 = percent_rank(CCL20),
         dist = rank_il6^2 + rank_ccl20^2) |>
  slice_max(dist, n = 11) |>
  mutate(type = ifelse(str_starts(name, 'Pept'), 'N', 'C'),
         tpp = str_extract(name, 'PTD4|R7|TAT'), id = str_extract(name, '\\d+$'),
         name = ifelse(name == '2-APB', '2-APB', str_c(tpp, '_', type, id)))

a431_cytk |>
  ggplot(aes(IL6, CCL20)) +
  geom_vline(xintercept = 1, linetype = 'dashed') +
  geom_hline(yintercept = 1, linetype = 'dashed') +
  geom_point() +
  geom_point(data = high_cytk, color = 'red') +
  geom_text_repel(data = high_cytk, aes(label = id)) +
  expand_limits(x = 0, y = 0) +
  labs(title = 'TRPV3-binding peptides stimulation on A431',
       x = 'Fold change of IL6', y = 'Fold change of CCL20') +
  theme_bw()

# A431 calcium flow ------------
cafl <- read_csv('virtual_screen_v3/data/a431-ca3.csv')

cafl_tidy <- cafl |>
  pivot_longer(-1) |>
  mutate(conc = str_remove(name, '.+-'),
         id = str_extract(name, '.+(?=-)'), name = NULL)

cafl_anno <- cafl_tidy |>
  left_join(a431_cytk) |>
  mutate(type = ifelse(str_starts(name, 'Pept'), 'N', 'C'),
         tpp = str_extract(name, 'PTD4|R7|TAT'), num = str_extract(name, '\\d+$'),
         name = ifelse(is.na(name), id, str_c(tpp, '_', type, num)),
         group = str_glue('{name} ({conc})') |> fct_reorder(value, .fun = max, .desc = T))

cafl_anno |>
  ggplot(aes(time, value, color = group)) +
  geom_path() +
  theme_pubr(legend = 'right') +
  labs(x = 'Time (s)', y = 'Fluo-4 AM', title = 'Ca2+ in A431', color = '') +
  scale_color_manual(values = c('red2', 'gold', 'darkgreen', 'royalblue',
                                grey.colors(9)))

cafl_peak <- cafl_anno |>
  slice_max(value, n = 3, by = group)

cafl_peak |>
  mutate(group = fct_reorder(group, -value)) |>
  ggplot(aes(group, value, fill = group)) +
  stat_mean(geom = 'col') +
  stat_summary(geom = 'errorbar', fun.data = mean_se, width = .5) +
  theme_pubr(legend = 'right') +
  labs(x = 'Group', y = 'Fluo-4 AM', title = 'Ca2+ peak in A431', fill = '') +
  scale_fill_manual(values = c('red2', 'gold', 'darkgreen', 'royalblue', grey.colors(9))) +
  theme(axis.text.x = element_blank())


